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WIENER CHAOS VS STOCHASTIC COLLOCATION METHODS FOR LINEAR 
ADVECTION-DIFFUSION-REACTION EQUATIONS WITH MULTIPLICATIVE 

WHITE NOISE t 

ZHONGQIANG ZHANG, MICHAEL V. TRETYAKOV, BORIS ROZOVSKII, AND GEORGE E. KARNIADAKIS 


Abstract. We compare Wiener chaos and stochastic collocation methods for linear advection-reaction- 
diffusion equations with multiplicative white noise. Both methods are constructed based on a recursive 
multi-stage algorithm for long-time integration. We derive error estimates for both methods and com¬ 
pare their numerical performance. Numerical results confirm that the recursive multi-stage stochastic 
collocation method is of order A (time step size) in the second-order moments while the recursive 
multi-stage Wiener chaos method is of order A N + A 2 (N is the order of Wiener chaos) for advection- 
diffusion-reaction equations with commutative noises, in agreement with the theoretical error estimates. 

However, for non-commutative noises, both methods are of order one in the second-order moments. 
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Notation. 

q: number of Brownian motions (noises). 

N: highest order of Hermite polynomial chaos. 

n: number of basis modes in approximating the Brownian motion. 

L: level of Smolyak sparse grid collocation. 

M: number of Fourier collocation nodes in physical space. 

A: element size (in time) for multi-element spectral approximation of Brownian motion. 

K: number of elements in time, which is T/A with T the final integration time. 

St: time step size for time discretization in the time interval (0, A]. 

? 7 (L, n q): number of sparse grid points at level L with dimension n q. 

1. Introduction 

Partial differential equations (PDEs) driven by white noise have different interpretations of stochas¬ 
tic products and lead to different numerical approximations, unlike the PDEs driven by colored noise. 
Specifically, stochastic products for white noise are usually interpreted with two different products: the 
Ito product and the Stratonovich product, see e.g. [1]. Though a problem can be equivalently formulated 
using these two products, the use of different products leads to different performance of numerical solvers 
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for PDEs driven by white noise, especially when Wiener chaos expansion (WCE) and stochastic colloca¬ 
tion methods (SCM) in random space are used. In this paper, we will show theoretically and through 
numerical examples that for white noise driven PDEs, WCE and SCM have quite different performance 
when the noises are commutative. This is different from how WCE and SCM behave for PDEs driven by 
colored noise. For elliptic equations with colored noise, it is demonstrated in [3, 10] that there are only 
small differences in the numerical performance of generalized polynomial chaos expansion and SCM. 

To apply WCE and SCM, we first discretize the Brownian motion with its truncated spectral expansion, 
see e.g. [29, Chapter IX] and [21], which results in PDEs with finite dimensional random inputs. Hence, 
our methods are Wong-Zakai type approximations [34, 35], where the Brownian motion is approximated by 
a smooth stochastic process of bounded variation, e.g., the spectral approximation used here and piecewise 
linear approximation of the Brownian motion [34, 35]. We note that piecewise linear approximation can 
be used instead but this is beyond the scope of the paper. 

The resulting PDEs can be solved numerically using a variety of space-time discretization methods 
and any sampling methods or functional expansion methods in random space. In random space, we will 
employ functional expansion methods, WCE [4, 21] and SCM [38], instead of the Monte Carlo method. 
These functional expansion methods have no statistical errors as no random number generators are used; 
they have only errors from truncations of Wiener processes and functional expansions and allow efficient 
short-time integration of SPDEs [4, 5, 15, 21, 22, 37, 38]. 

In principle, we can employ any functional expansion, however, different expansions are preferred 
for different stochastic products because of computational efficiency. In practice, WCE is associated 
with the Ito-Wick product, see (2.9), as the product is defined with Wiener chaos modes yielding a 
weakly coupled system (lower-triangular system) of PDEs for linear equations. On the other hand, 
SCM is associated with the Stratonovich product, see (2.15), yielding a decoupled system of PDEs. 
These different formulations lead to different numerical performance as we demonstrate in Section 4; in 
particular, WCE can be of second-order convergence in time while SCM is only of first-order in time 
in the second-order moments for commutative noises. Further, when the noises serve as the advection 
coefficients, SCM can be more accurate than WCE when both methods are of first order convergence as 
the SCM (Stratonovich formulation) can lead to smaller diffusion coefficients than those for WCE (Ito 
formulation). 

However, a fundamental limitation of these expansion methods is the exponential growth of error 
with time and the increasing complexity as the number of random variables is increasing, generated 
by the discretization of the Brownian motion. To deal with this complexity, a recursive WCE method 
was proposed in [21] for the Zakai equation of nonlinear filtering with uncorrelated observations. More 
recently, a recursive multi-stage approach was developed to efficiently solve linear stochastic advection- 
diffusion-reaction equations using WCE [37] or SCM [38]. 

To deal with the complexity in random space, some preprocessing procedures have been proposed, 
see e.g. [8, 31]. In these procedures, we are seeking the solution in the form u(t,x;u)) = E[u(t, x; •)] + 
i»(ij w)iq(f, x). Then by imposing the spatial orthogonality of Ui(x, t) and dt‘Uj{t, x) (i, j = 1,2,...), 
we can obtain an equivalent systems of SPDEs: a PDE for E[it(f, x\ •)], a system of equations for Yi{t,u>) 
and a system of equations for m(t , x). In many applications, this procedure is efficient, even with few terms 
of Yi(t.uj) and Ui(t,x ) as it may take advantage of some intrinsic sparsity structures of the underlying 
problems. However, this procedure requires some numerical methods to obtain Yi(t,u>), such as WCE 
(e.g. in [8]) and Monte Carlo methods (e.g. in [31]). When WCE or SCM is used, the complexity 
in random space is still high and thus the procedure is not efficient for problems driven by Brownian 
motion, where many modes Yi(t,u>) and Ui(t,x) are required. Though these procedures can be applied 
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here, we limit ourselves to the issue of using the deterministic integration methods, without using these 
procedures. 

Some numerical results of WCE for SPDEs have been presented in [37] for linear advection-diffusion- 
reaction equations and in [15] for nonlinear SPDEs including the stochastic Burgers equation and the 
Navier-Stokes equations. Numerical results for SCM have also been provided in [38] for linear stochas¬ 
tic advection-diffusion-reaction equations and the stochastic Burgers equation. Numerical results have 
demonstrated that WCE [37] and SCM [38] in conjunction with the recursive multi-stage approach are 
efficient for long-time integration of linear advection-diffusion-reaction equations. 

The main aim of the current paper is the derivation of theoretical error estimates for both WCE 
and SCM methods and subsequent comparison of the numerical performance of the two methods for 
commutative and non-commutative noises. In addition, we will develop a recursive multi-stage SCM, 
different than in [38], using a spectral truncation of Brownian motion. Specifically, in this paper we 
will derive the error estimate of WCE for linear advection-diffusion-reaction equations with white noise 
in the advection velocity and that of SCM with white noise in the reaction rate. We note that the 
convergence rate of WCE is known only for linear advection-diffusion-reaction equations with white noise 
in the reaction rate although the convergence of WCE for linear advection-diffusion-reaction equations 
has been studied for some time [20, 21, 22, 23]. 

The paper is organized as follows. After the Introduction section, in Section 2, we review the WCE 
and SCM for linear parabolic SPDEs and develop a new recursive SCM using a spectral truncation of 
Brownian motion, following the same recursive procedure as WCE in [21, 37, 38]. In Section 3 we present 
the error estimates for both methods for linear advection-diffusion-reaction equations, with the proofs 
presented in Section 5. In Section 4, we present numerical results of WCE and SCM for linear SPDEs 
with both commutative and non-commutative noises and verify the error estimates of WCE and SCM 
for commutative noises. 

2. Review of Wiener chaos and stochastic collocation 

In this section, we briefly review WCE and SCM for the following linear SPDE in the Ito form: 

9 

du{t,x) = £u(t,x)dt + y^Mku(t,x) dwk(t), (t,x) e (0,T] xD, 

k= 1 

u(0,x) = Uq(x), x€T>, ( 2 - 1 ) 

where = ({u’fc(f), 1 < k < q} ,Tt) is a system of one-dimensional independent standard Wiener 

processes defined on a complete probability space (fl, T, P) and 

d d 

Cu(t,x) = y aj t j (x) DjDju(t , x) + ^ bj(x)Dju(t, x) + c ( x ) u(t, x), 
i,j= 1 2—1 

d 

Mku(t,x) = i,k(x)Diu(t,x) + v k (x)u(t,x), (2.2) 

2=1 

and Di is the spatial derivative in a^-direction. We assume that the domain T> in R d is such that the 
periodic boundary conditions can be imposed or that T> = M d . In the former case, we will consider periodic 
boundary conditions and in the latter the Cauchy problem. 

We assume that there exist a constant Sc > 0 and a real number Cc such that for any v £ H 1 (T>), 

{Cv, v) + i ll-^HI 2 + 5 C IMIffi < Cc IMI 2 , 

k= 1 


(2.3) 
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where (•,•) is the duality between the Sobolev spaces i7 _1 (P) and H 1 (T>) associated with the inner- 
product over L 2 (V) and ||-|| is the P 2 (P)-norm. Specifically, we require that the coefficients of operators 
C and A4 in (2.2) are uniformly bounded and that 

d ( q 

2 a itj (x) - ^ cr z,k{ x ) a k,j( x ) 

i,j =1 \ k—1 

in addition to the Lipschitz continuity of a,ij(x). If E[||uo|| 2 ] is bounded (E[-] is the expectation with 
respect to P), these assumptions are sufficient for a unique square-integrablc solution of (2.1)-(2.2), see 
e.g. [23, 25], 

The problem (2.1)-(2.2) is said to have commutative noises if 



^ ViVj > \yf , x,y£V, 


M k Mj = MjM k , 1 <k,j<q, (2.4) 

and to have non-commutative noises otherwise. When q = 1, (2.4) is satisfied and thus this is a special 
case of commutative noises. When Mk are zeroth-order operators, (cr,;.*, = 0), (2.4) is satisfied and the 
problem also has commutative noises. The definition is consistent with that of commutative and non- 
commutative noises for stochastic ordinary differential equations, see e.g. [26]. In Section 4, we test our 
algorithms on examples with both commutative and non-commutative noises. 


Remark 2.1. The problem (2.1)-(2.2) can be regarded as an approximation of a problem driven by a 
cylindrical Wiener process. Consider a cylindrical Wiener process W(t,x) = ^kWk(t)ek{x) where 

EfcLi < oo, {uik(t)} are independent Wiener processes, and {efe(x)}^l 1 is a complete orthonormal 
basis (CONS) in L 2 (V), see e.g. [9, 30]. Thus, we can view (2.1)-(2.2) as approximations of SPDEs 
driven by this cylindrical Wiener process. 


In both WCE and SCM, we discretize the Brownian motion using the following spectral representations 
(see e.g. [21, 37]): 

n 

lim E[(w{t) - w^(t)) 2 } = 0, w {n \t) = Y & 

n—>-oo ^ J 

2=1 

where are mutually independent standard Gaussian random variables and is a CONS in 

P 2 ([0,P]). The expansion (2.5) is an extension of Fourier expansion of Brownian motion that is the 
Wiener construction [29, Chapter IX], see also [17, 18]. 


TOj(s) ds, t £ [0, T], 


(2.5) 


2.1. Wiener chaos expansion (WCE). The WCE solution to (2.1) is defined with the Cameron- 
Martin basis [6] in Wiener chaos space, using Fourier-Hermite series. The corresponding coefficients are 
obtained by solving the associated propagator , which is a lower-triangular linear system of deterministic 
parabolic equations determined by (2.1). Specifically, the solution to (2.1) can be represented as 


(t,x) = ^2 -/=^Ta(t,x-,u 0 )^ a , t£(0,T], 


a£j q 


( 2 . 6 ) 


where J q is the set of multi-indices a = (ak,i)k,i >l of finite length, i.e., 


Jq = < a = ( a k ,i , 1 < k <q), l > 1, a k ,i G {0,1,2,...}, |a| := a k ,i < oo > . 


k,l 


The random variables £ a are Cameron-Martin orthonormal basis, defined as 


n 


Hak'i (£fc,z) 


, a £ fj , 


(2.7) 
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where £^2 = / mi(s) dwk(s), and H n is the n-th Hermite polynomial: 

Jo 

H n { x) = (-0 n e x2/2 £^e- x212 . (2.8) 

Under our assumptions, the SPDE (2.1) can be written in the following form using the Ito-Wick product, 
see e.g. [14, Section 2.5] and [22], 

9 

du(t, x) = £u(f, x) dt -f ^2 Mku(t, x) o dw k {t ), (t, x) e (0, T] x D, 

k =1 

u(0,x) = uq(x), x £ V, (2-9) 


where the Ito-Wick product “o” is defined for the Cameron-Martin basis (2.7) such that o = 


a + /3)! 
al/31 


£a+/ 3 - By (2.6) and Cameron-Martin theorem [6], we obtain the coefficients <p a (t,x;uo) = 


E[v / alu(t,x)^ ct ]. Approximating w k with in (2.5), we substitute the representation (2.6) into (2.9) 
and then we can readily check that the coefficients ip a (t,x;) from (2.6) satisfy the following propagator, 
see e.g. [20, 22], 


dip a (t,x;u 0 ) 

dt 

<p a {0,x) 


q n 

Cip a (t,x;u 0 ) + EE a k ,imi(s)M k iPa-(k,i)(x;u 0 ), si E (0,T], 

k= 1 1=1 

M o(*)l{|a|=0}j 


where a (fc, l) is the multi-index with components 


(a (k,l)) 


max(0, aij — 1), if i = k and j = l, 
aij , otherwise. 


( 2 . 10 ) 


In practical computations, we have to truncate the propagator (2.10) and, consequently, we are inter¬ 
ested in the following truncated Wiener chaos solution: 


u N , n (t,x) 


E ■ 

Ot(=.SFu,n,q 


'a! 


-(fi a (t,x;u 0 )£, a 


( 2 . 11 ) 


where the set Jn,n, q = {a = {a k ,i)qx n | EaUi E"=i — N}. Here N is the highest Hermite polynomial 
order and n is the maximum number of Gaussian random variables for each Wiener process. In (2.7), we 
choose the basis {mz(s)} ;>1 as 


mi (s) = -E, mi(s) = cos(ELE£), / > 2, 0 <s<T. (2.12) 

As shown in [4, 21, 37], the error induced by the truncation of Wiener chaos expansion grows exponen¬ 
tially with time and thus WCE is not efficient for long-time integration. To control the error behavior, 
[37] proposes a recursive WCE (see Algorithm 2.2 below) for computing the second moments, E[it 2 (t, a;)], 
of the solution of the SPDE (2.1). Specifically, we discretize the Brownian motion using the following 
spectral representation in a multi-element version, i.e., using K multi-elements [21, 37]: 

K n rt k At 

w ( A ’ n \t) = ^2^2 m^k{s)ds^k, t G [0,T], (2.13) 

fc= 1 »= 1 Jtk-lAt 

where 0 = t 0 < ti < • • • < t[< = T, t*, A t is the minimum of t*, = fcA and t, {rrii t k}°°^ 1 is a CONS 
in L 2 ([tfc,tfc+i]), and G : A- are mutually independent standard Gaussian random variables. After the 
truncation of Brownian motion, we can have a similar propagator as (2.10). Noticing the linear property 
and Markovian properties of the solution to (2.1), we take the solution at U_i as initial condition to 
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Figure 1. Illustration of the idea of recursive multi-stage approach for long-time inte¬ 
gration in [37]. 


QSt ■■■ A 2A . iA 

Inn . -uiL _i_i_ 


T = KA 


solve the solution over (tfc_i,tfc]- Thus, we can recursively compute the covariance matrix at with the 
covariance matrix at the time instant tfc_i. We then have the following algorithm for the second moments 
of the approximate solution; see Figure 1 for an illustration and [37] for the derivation. 

Algorithm 2.2 (Recursive multi-stage Wiener chaos expansion, [37, Algorithm 2]). Choose the algo¬ 
rithm’s parameters: a CONS {e m (x)} m >i and its truncation {e m (a;)}^[ =1 ; a time step A; N and n which 
together with the number of noises q determine the size of the multi-index set Jn. n , q . 

Step 1. For each m = 1,..., M, solve the propagator (2.10) for a G fFi.n.q on the time interval [0, A] 
with the initial condition <f(x) = e m (x) and denote the obtained solution as ip a (A, x; e m ), rn = 1,...,M, 
and a G Jb.n.g Also, choose a time step size St to solve numerically the equations in the propagator. 

Step 2. Evaluate q a ,i,m = (y>a(A, •; e/), e m (-)), l,m = 1,..., M. Here (•, •) is the inner product in L 2 (V). 
Step 3. Recursively compute the covariance matrices Qi m (Ui N, n, M), l,m = 1,...,M, tj = iA, as 
follows: 

bl, M) — (^0: ei)(ll o, 6m); 

M 1 

Qlmifii N, fl, M) — ^ ] Q N, C. M) ^ f qa,j,iqa,k,mi 'l — 1) • * ■ ) K, 

j,k= 1 

where uq(x) is the initial condition for (2.1) and obtain N n (tj,x), the second moments of the approx¬ 
imate solution to (2.1) by the following: 

M 

M A,N,n(ti,a:) = Qim(ti;N,n,M)ei{x)e m (x), i = 1,..., K. (2.14) 

l,m—l 

Remark 2.3. The complexity of this algorithm is of order M 4 but can be reduced to the order of M 2 by 
making full use of the sparsity of the data [37]. 

2.2. Stochastic collocation method (SCM). This method leads to a fully decoupled system instead 
of a weakly coupled system from the WCE. First, we rewrite the SPDE (2.1) in the Stratonovich form, 
see e.g. [13, 19], 

9 

du(t, x) = Cu{t, x) dt + ^ A iku(t,x) o dwk{t), {t,x) G (0,T] x V, 

fc=l 

u(0,a;) = uo(x), xGV, (2-15) 

where Cu = Cu — |X^i<fc< 9 AA k A4kU. Second, we approximate the Brownian motion with its multi¬ 
element spectral expansion (2.13), and obtain the following partial differential equation with smooth 
random inputs (see e.g. [13]): 

9 

dUA.n(t,x) = CuA.n(t,x) dt + y^ y MkUA, n {t,x)dw ( ' k A ’"\t), (t,x) G (0, T] X V , 

fc=1 


u(0, x ) 


uq{x), x G V . 


(2.16) 
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In (2.16), we have ngK standard Gaussian random variables j, l < n, k < q, i < K, according to (2.13). 
Now we can apply standard numerical techniques of ngK-dimensional integration to numerically obtain 
p-th moments of the solution to (2.16): 

E[(uA,n(7>)) p ] = —— l —^ [ (F(u 0 (x),T,x,y)) p e-^ dy, p = 1,2,... (2.17) 

where y = ( yi t k,i ), l < n, k < q, i < K and the functional F represents the solution functional for (2.16). 
Here we use sparse grid collocation [12, 33] if the dimension ngK is moderately large. As pointed out in 
[2, 36], we are led to a fully decoupled system of equations as in the case of Monte Carlo methods. 

In practice, we use the following sparse grid quadrature rule for a d-dimensional function ip, see e.g. 
[12, 33], 

A(L,d)p= Y, (-l) L+d - 1 - |i ^“MQ il ®--.®Q irf p, (2.18) 

L<|i|<L+d-l ^ ' ’ 

where we have one-dinrensional Gauss-Hermite quadrature rules Q n for univariate functions i/>(y), y £=R: 
QnV'(y) = )Ca=i ip{yn,a)wn,a, Yn,i < Yn, 2 < < y n ,n are the roots of the n-th Hermite polynomial 

(2.8) and w „ |0 are the associated weights w njQ = n!/n 2 [77 n _i(y„ iQ ,)] 2 . The number of sparse grid points, 
denoted by t/(L, d), for this sparse grid quadrature rule is of order d L_1 when L < d, which can be checked 
readily from the rule (2.18). For example, we have, for L = 2,3,4, 

77 ( 2 , d) = 2d + 1, p(3, d) = 2d 2 + 2d + 1, 77 ( 4 , d) = -d 3 + 2d 2 + ^d + 1. 

o o 

Denote the set of ? 7 (L, d) sparse grid points x„, = (x)., • • • , x(?) by H]] 9 , where (1 < j < d) belongs to the 
set of points used by the quadrature rule . According to (2.18), we only need to know the function 
values at the sparse grid H 

’j(L.nq) 

A{L,d)p= Y ^( X «) W «J x « = > x «) € n n L q , (2.19) 

K,= 1 

where W K are determined by (2.18) and the choice of the quadrature rules Qi j and they are called the 
sparse grid quadrature weights. 

Here again, the direct application of SCM is efficient only for short-time integration. To achieve 
long-time integration, we apply the recursive multi-stage idea used in Algorithm 2.2, i.e., we use SCM 
over small time interval (tj_i,tj] instead of over the whole interval (0,T] and compute the second-order 
moments of the solution recursively in time. The derivation of such a recursive algorithm will make use 
of properties of the problem ( 2 . 1 ) and orthogonality of the basis both in physical space and in random 
space as will be shown shortly. 

We solve (2.16) with spectral methods in physical space, i.e., using a truncation of a CONS in physical 
space {e m }^ =1 to represent the numerical solution. The corresponding approximation of u& tn (t,x) is 
denoted by u^ n (t,x). Further, let v(t,X’ : s,vo) be the approximation u^ n (t,x) of UA,n(t,x) with the 
initial data vq prescribed at s : UA,n(s,a:) = vq{x). Note that 

MA,n(ti, z) = v{t U X- ti-l, WA,n(ti-l, '))> U = lA. (2.20) 

Denote 4> m (ti; A, n, M) = (u^ n (tj, •), e m ). Then the second moments are computed by 

M 

E [(“A,n( t i) a: )) 2 ] = Y H im(U;A,n,M)ei(x)e m (x), 

Z,ra—1 


( 2 . 21 ) 
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where A, n, M) = E[$j(tj; A, n, M)<f> m (ti; A, n, M)]. Now we show how the matrix Hi m (tp 7 A, n, M) 

can be computed recursively. By the linearity of (2.16), we have 


M 


WA.nftti 1 ) = Y. ®i{ t»-i; A, n, M)u(h, x; h_i ,ej). 


Denote " ( V {U, _i, ef), e m ). Then by the orthonormality of e m , we have 


M 



The matrix H[ rn (t l ; A, n, M) can be computed recursively as 


M M 



j =1 fc=l 


We note that the expectation F.[hj : i t i-ihk,m,i-i] does not depend oni-1 because according to (2.16) and 
(2.2), v(U,x-,U-i ,ei) depend on the length of the time interval A and the random variables (l < n, 
k < q) but is independent of time t*_i- Denote v(ti, ■; tj_i, e{) with £i t k,i anchored at the sparse grid 
point x K £ Vff by v K (A,-\ei). Let h K ,i,m = (v K (A, •; ei), e m ). Then, using the sparse grid quadrature 
rule (2.19), we obtain the recursive approximation of if; m (tj; A, n, M): 


»7(L,nq) 


M M 



Him(U\ A, n, M) « A, L, n, M) := ^ ^ A, L, n, M) 


j = l k =1 


Substituting (2.22) in (2.21), we obtain an approximation for the second moments of u(t 7 x) 7 denoted by 
L n (t i , x ). When M = oo (i.e., when the CONS { e m } is not cut-off), we denote this approximation by 

M A ,L.n(ti, X). 

Remark 2.4. For non-homogeneous equations, i.e., with forcing terms, we can have similar algorithms. 
Indeed, the same procedure applies once we can split the non-homogeneous equations into two equations: 
non-homogeneous equation with zero initial value and homogeneous equation with initial value. See [38] 
for a derivation of similar algorithms where only increments of Brownian motion are used, which is 
different from the spectral approximation of Brownian motion used here. 

Now we have the following algorithm for the second moments of the approximate solution. 

Algorithm 2.5 (Recursive multi-stage stochastic collocation method). Choose a CONS {e m (x)} m >i 
and its truncation {e m (:r)}^ =1 ; a time step A; the sparse grid level L and n, which together with the 
number of noises q determine the sparse grid 'H'f 1 which contains r/(L. n q) sparse grid points. 

Step 1. For each m = 1,.... M, solve the system of equations (2.16) on the sparse grid "H]] 9 in the time 
interval [0, A] with the initial condition <j>(x) = e m (x) and denote the obtained solution as u K (A,x;e m ), 
to = 1,.... M, and k = 1,- ■ ■ 7 r]( L, ng). Also, choose a time step size 5t to solve (2.16) numerically. 

Step 2. Evaluate h K j t m = (w K (A, -;ej), e m ), l,m = 1,..., M. 

Step 3. Recursively compute the covariance matrices Hi m ( tjj L, n, M), Z,m=l,...,M, as follows: 


A, L, n, M) (lip , e^) (l/Q; ^-m) 


M 


»j(L,n (?) 


Him(t i7 A, L, n,M) 
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where Ug(x) is the initial condition for (2.1) and obtain the approximate second moments L n (t i,x) 

of the solution u(t,x) to (2.1) as 

M 

= ^ Hj m (ti; A, L, n, M)e ; (a;)e m (a;), i = 1,..., K. (2.23) 

l,m=l 

Remark 2.6. Similar to Algorithm 2.2, the cost of this algorithm is ^? 7 (L.ng)M 4 and the storage is 
r)(L, ng)M 2 . The total cost can be reduced to the order of M 2 by adopting reduced order methods in 
physical space, see e.g. [32]. The discussion on computational efficiency of the recursive WCE methods, 
see [37, Remark 4.1], is also valid for Algorithm 2.5. 


3. Error estimates 


Though WCE and SCM use the same spectral truncation of Brownian motion, the former is associated 
with the Ito-Wick product while the latter is related to the Stratonovich product. Note that WCE employs 
orthogonal polynomials as a basis and SCM does not have such orthogonality. This difference allows WCE 
to have a better convergence rate than SCM in the second-order moments, see Corollary 3.2 and Theorem 
3.3. 

Assume that the operator C generates a semi-group {7t} t>0 , which has the following properties: for 
v G H k (V), 

||7tv||^ fc < C(k, C)e 2Cct , (3.1) 

where C(0,£) = 1 and 

[ e 2C ^-^ \\T t v\\ 2 Hk+1 dO < S^C{k, C)e 2Gc< ' t ~ s ' > \\vf Hk . (3.2) 

J S 

Also, we assume that there exists a constant C[r,M.) such that 

\\Mivf Hk < C(k,M) \\v\\ 2 Hk+ i , for v G H k+1 , l = l,...,q. (3.3) 

and that there exists a constant C(k , C) such that 

\\Cv\\ 2 Hk < C{k,C) \\v\\ 2 Hk+2 , for v G H k+2 . (3.4) 


The conditions (3.1) and (3.3) are satisfied with k < r and (3.4) is satisfied with k < r — 1 when the 
coefficients from (2.2) belong to the Holder space Cf +1 (V), which is equipped with the following norm 


max \\Dnf\\ Ta 
0<|/3|<LrJ 11 PJ " L 


sup 

x,y£T> 

|/9|=LrJ,r>|rJ 


\Dpf{x) - Dpf(y )| 


' - y I 


r— [r*J 


and [tJ is the integer part of the positive number r, cf. [11, Section 5.1]. Define also 

C k = max |c(j, C)C(j - 1, A4)} . 


(3.5) 


For the WCE for the SPDE (2.1) with single noise (q = 1), we have the convergence results stated 
below. In the general case, we have not succeeded in proving such theorems but we numerically check 
convergence orders using examples with commutative noises and non-commutative noises in Section 4. 


Theorem 3.1. Let q = 1 in (2.1). Assume that in (2.2) belong to Cf +1 (T>) and uo G 

H r {T>), where r > N + 2 and N is the order of Wiener chaos. Also assume that (2.3) holds. Then for 
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C i < 6c, the error of the truncated Wiener chaos solution Un n (tj, x) from (2.11) is estimated as 

( E [|| UN , n ( t i ,-)-«( t ij -)|| 2 ]) 1/2 

(C7 Lr jA)W-N-i Sc 


e C M T 

(N + i)! ^ LrJT 

A IKII 


< (C LT .jA) N / 2 e CcT 


-^2C n+ 2 C(N +2,C)C{NX)e Cfl+2T+CcT 


1/2 


IKII 


H r 


\/r\7T 


£[H+2 , 


(3.6) 


where t* = iA, the constants Sc and Cc are from (2.3), C i r j is defined in (3.5), £7(1X1,/!) is from (3.4), 
and £7(N + 2, C) is from (3.1). 


From Theorem 3.1, we have that the mean-square error of the recursive multi-stage WCE is 0(A N / 2 ) + 
0(A). The same result is proved for q = 1 and crj >r = 0 in [21], where the condition C\ < Sc is not 
required. Also, for the case of a^ r 0, the mean-square convergence without order but not requiring the 
condition C\ < Sc was proved in [20, 22]. 


Corollary 3.2. Under the assumptions of Theorem 3.1, we have 

-E[|K ti ,.)|| 2 ]|=E[||u N , 
e c M T (C LrJ A)W-N-i 


E[||uN,n(ti, -)l| 2 ] - EflKti, -)|| 2 ] = E[||uN, n (tj, •) - u(u, -)|| 2 ] 


< (q,j A ) e 


N2CcT 


(N + l)! 


LrJ! 


Sc - C\ 


IKII*, 


+2£7 n+2 £7(N +2,C)C(N,C)e 2CN+2T+2CcT ^ |KII*~+» . 

n7r 


(3.7) 


This corollary states that the convergence rate of the error in second-order moments (3.7) is twice that 
of the mean-square error (3.6) , i.e., 0(A N ) + 0(A 2 ). This corollary can be proved by the orthogonality 
of WCE. In fact, it holds that 


E[u 2 (ti, a;)] - E[i4, n (ti, x)] = E[(u(t,, x) - u N .n(ti, x)) 2 ], (3.8) 

as the different terms in the Cameron-Martin basis are mutually orthogonal [6]. Then integration over 
the physical domain and by the Fubini Theorem, we reach the conclusion in Theorem 3.1. 

For SCM for the SPDE (2.1), we have the following estimates: the first one is weak convergence of 
the Wong-Zakai type approximation UA, n (t,x) from (2.16) to u(t,x ) from (2.1), see Theorem 3.3; the 
second one is the convergence of SCM, i.e., the convergence of MA,L.n(t», x) to E[u^ n (t», x)], see Theorem 
3.4. Here we prove the convergence rate when Oi >r = 0 which belongs to the case of commutative noises 
(2.4). Our proof for Theorem 3.3 is based on the mean-square of convergence of the Wong-Zakai type 
approximation (2.16) to (2.1). When cq jr ^ 0, we have not succeeded in proving this mean-square 
convergence and, as far as we know, only a rate of almost sure convergence of the Wong-Zakai type 
approximations to (2.1) has been proved so far [13]. 


Theorem 3.3. Assume that Gi^ r = 0 and that the initial condition uq and the coefficients in (2.2) are in 
Cfr(D). Let u(t,x) be the solution to (2.1) and UA, n (t,x) be the solution to (2.16). Then for any e > 0, 
there exists a constant C > 0 such that the one-step error is estimated by 

|E[u 2 (A,x)] -E[u^ n (A,x)]| < £7exp(CA)(A 6 +A 2 )n- 1+£ , (3.9) 

and the global error is estimated by 

|E[u 2 (ti,x)] -E[{& in (ti,x)]| < £7exp(£7T)ArW 1+E , 1 < * < K. (3.10) 

The following theorem is on the convergence of the second moments by SCM to those of the solution 
to (2.16). 
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Theorem 3.4. Let u^ n (t,x) be the solution to (2.16) and M Ai |_. n (tj, x)be the limit of M“ l n (tj, x) from 
(2.23) when M — > oo. Under the assumptions of Theorem 3.3, for any e > 0, the one-step error is 
estimated by 

|M A , L . n (A,x) - E[ui in (A,z)]| < C exp(CA)(A 3L + A 2L )(1 + (3c/2) LAn )/r0- An )/ 2 £ - L L- 1 n u , 

and the global error is estimated by, for 1 < * < K, 

|M AjL)I1 (U, x') E[rt A „(t A x)]| <Cexp(CT)A 2L - 1 (l + (3c/2) LAn )/3-( LAn )/ 2 £- L L- 1 n u . 

Here the positive constants C, c, /3 < 1 are independent of A, L and n. The expression L A n means the 
minimum of L and n. 

According to Theorems 3.3 and 3.4, the error of the SCM is 0(A 2L_1 ) + 0(A) in the second-order 
moments. Compared to Corollary 3.2, the SCM is of one -order lower than WCE when N = 2 as the 
error of WCE is 0(A N ) + 0(A 2 ). 


4. Numerical results 


In this section, we compare Algorithms 2.2 and 2.5 for linear stochastic advection-diffusion-reaction 
equations with commutative and non-commutative noises. We will test the computational performance 
of these two methods in terms of accuracy and computational cost. All the tests were run using Matlab 
R2012b, on a Macintosh desktop computer with Intel Xeon CPU E5462 (quad-core, 2.80 GHz). Every 
effort was made to program and execute the different algorithms as much as possible in an identical way. 

We note that we do not have exact solutions for all examples and hence evaluate the errors of the 
second-order moments using reference solutions, denoted by E[u 2 ef (T, x)], which are obtained by either 
Algorithm 2.2 or Algorithm 2.5 with fine resolution. We do not use solutions obtained from Monte Carlo 
methods as reference solutions since Monte Carlo methods are of low accuracy and are less accurate than 
the recursive multistage WCE, see [37] for comparison between WCE and Monte Carlo methods and also 
below. 

The following error measures are used in the numerical examples below: 


e 2 (T) = |||E[ U 2 ef (T,.)]|| ; 2 -||M^(T,.)| 


^(T)=|||EK 2 ef (T,.)]|| ioo -||M^(T,.))| 


i 2 1 > 


Q r f\T) = 


el{T) 


l|Ek 2 ef (T, 


r ,OO /rji\ 

Qi ( T ) = 


Q?(T) 


Pkefk,-) 


(4.1) 


(4.2) 


where M^(T, x) is either N n (T, x) from Algorithm 2.2 or L n (T, x) from Algorithm 2.5, \\v\\ l2 = 
m \ i 


m E v2 m 


IH|;oo = max \v(x m )\, x m are the Fourier collocation points. 
l<m<M 


The computational complexity for Algorithm 2.2 is ( N j J n9 )- A -M 4 (see [37]) and that for Algorithm 2.5 
is rj(L, ng)-^M 4 . The ratio of the computational cost of SCM over that of WCE is ry(L, n q)/ ( N jj n<? ). For 
example, when N = 1 and L = 2, the ratio is (1 + 2n<?)/(l + n q), which will be used in the three numerical 
examples. The complexity increases exponentially with n q and L, see e.g. [12], or N but increases linearly 
with ^. Hence, we only consider low values of L and N. 


Example 4.1 (Single noise). We consider a single noise in the Ito SPDE (2.1) over the domain (0, T\ x 
(0,27r): 

du = [(e+-o 2 )d^.u +Psin(x)d x u]dt +ad x udw(t), (4.3) 

or equivalently in the Stratonovich form 

du = [ed 2 u +/3 sin(x)d x u] dt + ad x u o dw(t), 


(4.4) 
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with periodic boundary conditions and non-random initial condition u(0,x) = cos(x), where w(t) is a 
standard scalar Wiener process, e > 0, (3, a are constants. 

In this example, we compare Algorithms 2.2 and 2.5 for (4.3) with the parameters /3 = 0.1, a = 0.5 and 
e = 0.02. We will show that the recursive multi-stage WCE is at most of order A 2 in the second-order 
moments and the recursive multi-stage SCM is of order A. 

In Step 1, Algorithm 2.2, we employ the Crank-Nicolson scheme in time and Fourier collocation 
in physical space. We obtain the reference solution by Algorithm 2.2 with the same solver but finer 
resolution as a reference solution 1 since we have no exact solution to (4.3). The reference solution is 
obtained by M = 30, A = 10~ 4 , N = 4, n = 4, St = 10 -5 . It gives the second-order moments in Z 2 -norm 
||E[M 2 ef ]|| i2 =1.065194550063 and in the Z°°-norm ||E[u 2 ef ] || ;ao =0.5174746141105. 

From Table 1, we observe that the recursive WCE is 0(A N ) + 0(A 2 ) for the second-order moments. 
When N = 2, the method is of second-order convergence in A and of first-order convergence when N = 1. 
When N = 3, the method is still second-order in A (not presented here). This verifies the estimate in 
Corollary 3.2. 

Table 1. Algorithm 2.2: recursive multi-stage Wiener chaos method for (4.3) at T = 5: 


0.5,0 = 

0.1, e = 

0.02, and M = 

20, n = 

1. 



A 

St 

N 

e r i\T) 

order 

r.OO /rji\ 

02 [ T ) 

order 

CPU time (sec.) 

1.0e-l 

1.0e-2 

1 

1.5249e-2 


8.8177e-3 


3.57 

1.0e-2 

1.0e-3 

1 

1.5865e-3 

^0.98 

8.9310e-4 

^0.99 

33.22 

1.0e-3 

1.0e-4 

1 

1.5934e-4 

A 1.00 

8.9429e-5 

A 1.00 

348.03 

1.0e-l 

1.0e-2 

2 

1.9070e-4 


4.1855e-5 


5.14 

1.0e-2 

1.0e-3 

2 

2.0088e-6 

A 1.98 

4.2889e-7 

A 1.99 

51.75 

1.0e-3 

1.0e-4 

2 

2.0386e-8 

A 1.99 

4.8703e-9 

A 1 ' 94 

490.04 


In Step 1, Algorithm 2.5, we use the Crank-Nicolson scheme in time and Fourier collocation method in 
physical space. The errors are also measured as in (4.1) and (4.2). The reference solution is obtained by 
Algorithm 2.2 as in the case of WCE. We observe in Table 2 that the convergence order for second-order 
moments is one in A even when the sparse grid level L is 2, 3 and 4 (the last is not presented here). The 
errors for L = 3 are more than half in magnitude smaller than those for L = 2 while the time cost for 
L = 3 is about 1.5 times of that for L = 2. 

In summary, from Tables 1 and 2, we observe that the recursive multi-stage WCE is 0(A N ) + 0( A 2 ) 
and the recursive multi-stage SCM is 0(A), as predicted by the error estimates in Section 3. While the 
SCM and the WCE are of the same order when N = 1 and L > 2, the former can be more accurate than 
the latter. In fact, when N = 1 and L = 2, the recursive multi-stage SCM error is almost two orders of 
magnitude smaller than the recursive multi-stage WCE while the computational cost for both is almost 
the same, as predicted (( N ^ n9 ) = ry(L, n q) = 2). The recursive multi-stage WCE with N = 2 is of order 
A 2 and its errors are almost two orders of magnitude smaller than those by the recursive multi-stage 
SCM (with level 2 or 3) for the second-order moments. 

In this example, the recursive multi-stage SCM outperforms the recursive multi-stage WCE with 
N = 1. The reason can be as follows. In SCM, we solve an advection-dominated equation rather than a 

1 For single noise, it is proved in Theorem 3.1 that the recursive multi-stage WCE is of second-order convergence in 
second-order moments. The second-order convergence is numerically verified in [37]. For this specific example, a Monte 
Carlo method with 10® sampling paths (which costs 27.6 hours) gives | = 1.06517 ±6.1 X 10 —4 and || = 

0.51746±6.1 X 10“ 4 , where the numbers after ‘±’ are the statistical errors with the 95% confidence interval. We use Fourier 
collocation in space with M = 20 and Crank-Nicolson in time with St = 10“ 3 for the Monte Carlo method. 
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Table 2. Algorithm 2.5: recursive multi-stage stochastic collocation method for (4.3) 


= 5: a 

= 0.5, 0 

— 

0.1, e = 0.02, 

and M 

= 20, n = 1. 



A 

St 

L 

e?{T) 

order 

r, OO /rri\ 

Qi (T) 

order 

CPU time (sec.) 

le-01 

le-02 

2 

3.4808e-04 


3.0383e-03 


3.71 

le-02 

le-03 

2 

3.4839e-05 

A i.oo 

3.0130e-04 

A i.oo 

33.88 

le-03 

le-04 

2 

3.4844e-06 

A i.oo 

3.0106e-05 

A i.oo 

325.06 

le-01 

le-02 

3 

1.6815e-04 


3.4829e-04 


5.16 

le-02 

le-03 

3 

1.6230e-05 

A 102 

3.2283e-05 

A 1.03 

50.59 

le-03 

le-04 

3 

1.6170e-06 

A i.oo 

3.2026e-06 

a i.oo 

486.08 


diffusion-dominated equation in WCE, as SCM is associated with the Stratonovich product which leads 
to the removal of the term \a 2 df.u in the resulting equation, see (4.4). The larger cr is, the more dominant 
the diffusion is. In fact, results for cr = 1 and cr = 0.1 (not presented here) show that when a = 1, the 
relative error of SCM with L = 2 is almost three orders of magnitude smaller than WCE with N = 1; 
when cr = 0.1, the relative error of SCM with L = 2 is only less than one order of magnitude smaller than 
WCE with N = 1. With the Crank-Nicolson scheme in time and Fourier collocation in physical space, 
we cannot achieve better accuracy for WCE with N = 1 and A no less than 0.0005 when M < 40. 


Example 4.2 (Commutative noises). We consider two commutative noises in the Ito SPDE (2.1) over 
the domain (0, T] x (0, 27r): 

du = [(e-\-^(r\cos 2 {x))d 2 u-\-{f3svn{x)— ^a\s\a.(2x))d x u]dt 

+cri cos(x)d x udwi(t) + (J 2 udw 2 (t ), (4-5) 

or equivalently in the Stratonovich form 

du = [ed 2 u + /3sm(x)d x u] dt + or cos(x)d x u o dw\{t) + cr 2 u o dw 2 (t), (4.6) 

with periodic boundary conditions and non-random initial condition w(0, x) = cos(x), where W 2 (t)) 

is a standard two-dimensional Wiener process, e > 0, /?, cr 1; cr 2 are constants. The problem has commu¬ 
tative noises, see (2.4). 


In this example, we take 01 = 0.5, cr 2 = 0.2, /? = 0.1, e = 0.02. We again observe first-order convergence 
for SCM and WCE with N = 1, and second-order convergence for WCE with N = 2 as in the last example 
with single noise. 

We choose the same space-time solver for the recursive multi-stage WCE and SCM as in the last 
example. We compute the errors as in (4.1) and (4.2). In Table 3, the reference second moments are 


H M 

U A=10- 4 ,N 




and 


M^_iq_ 4 n n (T, •) obtained by Algorithm 2.2 with = 10 5 and all the 


other truncation parameters are the same as stated in the table. In Table 4, the reference second moments 


are 


M^o-.L.nCO 


and 


Mi =1 o-‘, L ,„( T ) 


obtained by Algorithm 2.5 with St = 10 5 while all the 
other truncation parameters are the same as in the table. 

Here we do not compare the performance of Monte Carlo simulations with our algorithms as the main 
cost of Monte Carlo methods is to reduce the statistical errors. For the same parameters described above, 
when we used 10 6 Monte Carlo sampling paths, we could only reach the statistical error of 8.3 x ICC 4 , 
in 3.9 hours. To obtain an error of 1 x 10 -5 , seven thousands times more Monte Carlo sampling paths 
should be used, requiring three years of computational time and thus is not considered here. In the next 
example, we have similar situations and hence we will not consider Monte Carlo simulations. This also 
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demonstrates the computational efficiency of Algorithms 2.2 and 2.5 in comparison with Monte Carlo 
methods. 

For WCE, we observe in Table 3 convergence of order A N (N < 2) in the second-order moments: 
first-order convergence when N = 1, and second-order convergence when N = 2. Numerical results for 
N = 3 (not presented here) show that the convergence order is still two even though the accuracy is 
further improved when N increases from 2 to 3. This is consistent with our estimate 0(A N ) + 0( A 2 ) in 
Corollary 3.2. 

We also tested the case n = 2 which gives similar results and the same convergence order. 

Table 3. Algorithm 2.2: recursive multi-stage Wiener chaos expansion for commutative 
noises (4.5) at T = 1: <j\ = 0.5, oi = 0.2, /3 = 0.1, e = 0.02, and M = 30, n = 1. 


A 

St 

N 

srm 

order 

-V, OO /rri\ 

Qi (T) 

order 

CPU time (sec.) 

le-01 

le-02 

1 

1.6994e-03 


1.6548e-03 


3.19 

le-02 

le-03 

1 

1.7838e-04 

^0.98 

1.7172e-04 

^0.98 

32.74 

le-03 

le-04 

1 

1.6323e-05 

A 1 - 04 

1.5694e-05 

A 1-04 

329.15 

le-01 

le-02 

2 

4.0658e-05 


2.9568e-05 


6.53 

le-02 

le-03 

2 

4.4805e-07 

A 1.96 

3.3106e-07 

A 1.95 

65.89 

le-03 

le-04 

2 

4.4682e-09 

^2.00 

3.3484e-09 

^2.00 

657.55 


For SCM, we observe first-order convergence in A from Table 4 when L = 2,3. We note that further 
refinement in truncation parameters in random space, i.e., increasing L and/or n do not change the 
convergence order nor improve the accuracy. The case L = 3 actually leads to a bit worse accuracy, 
compared with the case L = 2. We tested the case L = 4 which leads to the same magnitudes of errors as 
L = 3. We also tested n = 2 and observed no improved accuracy for L = 2,3, 4. These numerical results 
are not presented here. 

For the two commutative noises, we conclude from this example that the recursive multi-stage WCE 
is of order A N + A 2 in the second-order moments and that the recursive multi-stage SCM is of order A in 
the second-order moments no matter what sparse grid level is used. The errors of recursive multi-stage 
SCM are one order of magnitude smaller than those of recursive multi-stage WCE with N = 1 while the 
time cost of SCM is about 1.6 times of that cost of WCE. For large magnitude of noises (cri = a 2 = 1, 
numerical results are not presented), we observed that the SCM with L = 2 and WCE with N = 1 have 
the same order-of-magnitude accuracy. In this example, the use of SCM with L = 2 for small magnitude 
of noises is competitive with the use of WCE with N = 1. 

Table 4. Algorithm 2.5: recursive multi-stage stochastic collocation method for com¬ 
mutative noises (4.5) at T = 1 : cri = 0.5, 02 = 0.2, /3 = 0.1, e = 0.02, and M = 30, 
n = 1. 


A 

5t 

L 

Qi\T) 

order 

-r.OO /rji\ 

Qi (T) 

order 

CPU time (sec.) 

le-01 

le-02 

2 

1.3624e-04 


1.2453e-03 


5.18 

le-02 

le-03 

2 

1.3064e-05 

A 1-02 

1.2009e-04 

A 1 - 02 

54.70 

le-03 

le-04 

2 

1.1837e-06 

A 1-04 

1.0889e-05 

A 1 - 04 

545.20 

le-01 

le-02 

3 

2.5946e-04 


2.0482e-04 


13.26 

le-02 

le-03 

3 

2.5437e-05 

A 101 

1.7897e-05 

A 1.06 

142.23 

le-03 

le-04 

3 

2.3102e-06 

A 1-04 

1.6062e-06 

A 1.05 

1420.24 


























COMPARISON OF WIENER CHAOS AND STOCHASTIC COLLOCATION 


15 


Example 4.3 (Non-commutative noises). We consider two non-commutative noises in the Ito SPDE 
(2.1) over the domain (0,T] x (0, 27r): 

du = [(e + ^cr\)d x u + f3sm(x)d x u + -<r| cos 2 (:r)u] dt 

+a\d x udwi{t) + (72 cos(x)udw 2 (t), (4.7) 

or equivalently in the Stratonovich form 

du = [ed x u + /3 sin(a;)9 a: w] dt + u\d x u o dw\{t) + 02 cos (x)u o dw 2 {t), (4.8) 

with periodic boundary conditions and non-random initial condition u(0, x) = cos(a;), where (wi(f),W 2 (t)) 
is a standard Wiener process, e > 0, (3, or, 02 ore constants. The problem has non-commutative noises 
as the coefficients do not satisfy (2.4). 

We take the same constants e > 0, /?, or, 02 as in the last example. We also take the same space-time 
solver as in the last example. In the current example, we observe only first-order convergence for SCM 
(level L = 2, 3,4) and WCE (N = 1, 2, 3) when n = 1, 2, see Table 5 for parts of the numerical results. 

The errors are computed as in the last example. The reference solutions are obtained by Algorithm 
2.2 for the recursive multi-stage WCE solutions and by Algorithm 2.5 for the recursive multi-stage SCM 
solutions, with A = 5 x ICC 4 and St = 5 x 10 -5 and all the other truncation parameters the same as 
stated in Tables 5 and 6. 

Table 5. Algorithm 2.2 (recursive multi-stage Wiener chaos expansion, left) and Algo¬ 
rithm 2.5 (recursive multi-stage stochastic collocation method, right) for (4.7) at T = 1: 

(T 1 = 0.5, (T 2 = 0.2, /3 = 0.1, e = 0.02, and M = 20, n = 1. The time step size St is A/10. 

The reported CPU time is in seconds. 


A 

N 

&\T) 

order 

time (sec.) 

L 

e7(T) 

order 

time (sec.) 

1.0e-l 

1 

3.7516e-03 


1.04 

2 

6.4343e-04 


1.65 

5.0e-2 

1 

1.8938e-03 

A o.99 

2.11 

2 

3.1738e-04 

A 102 

3.31 

2.0e-2 

1 

7.5292e-04 

A 1 - 01 

5.12 

2 

1.2440e-04 

A 102 

8.64 

1.0e-2 

1 

3.6796e-04 

A 1.03 

10.19 

2 

6.0502e-05 

A 104 

17.12 

5.0e-3 

1 

1.7457e-04 

A 1.08 

20.01 

2 

2.8635e-05 

A i.os 

33.82 

2.0e-3 

1 

5.8246e-05 

A 1 ' 20 

50.39 

2 

9.5401e-06 

A 1 - 12 

86.44 

1.0e-l 

2 

9.4415e-05 


2.16 

3 

1.5803e-04 


4.03 

5.0e-2 

2 

3.7303e-05 

A 1 - 81 

4.11 

3 

7.6548e-05 

A 1.05 

8.68 

2.0e-2 

2 

1.2282e-05 

A 1 - 34 

9.97 

3 

2.9673e-05 

A 1.03 

22.08 

1.0e-2 

2 

5.5807e-06 

A 1 - 21 

20.03 

3 

1.4378e-05 

A 1.05 

43.85 

5.0e-3 

2 

2.5471e-06 

A 1 - 14 

40.25 

3 

6.7925e-06 

A 1.08 

88.35 

2.0e-3 

2 

8.2965e-07 

A 1 - 22 

101.34 

3 

2.2605e-06 

A 1-20 

223.15 


In this example, our error estimate for recursive multi-stage WCE is not valid any more and the 
numerical results suggest that the errors behave as A N + C A/n. For N = 1 and n = 10 (not presented), 
the error is almost the same as n = 1. While N = 2 and n = 10, the error first decreases as 0{ A 2 ) for 
large time step size and then as O(A) for small time step size; see Table 6. When N = 2 and n = 10, the 
errors with A = 0.005,0.002,0.001 are ten percent (1/n) of those with the same parameters but n = 1 
in Table 5. Here the constant in front of A, C/n, plays an important role: when A is large and this 
constant is small, then the order of two can be observed; when A is small, C A/n is dominant so that 
only first-order convergence can be observed. 
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The recursive multi-stage SCM is of first-order convergence when L = 2,3,4 and n = 1, 2,10 (only parts 
of the results presented). In contrast to Example 4.2, the errors from L = 3 are one order of magnitude 
smaller those from L = 2. Recalling that the number of sparse grid points is rj( 2, 2) = 5 and r)(3, 2) = 13, 
we have the cost for L = 3 is about 2.6 times of that for L = 2. However, it is expected that in practice, 
a low level sparse grid is more efficient than a high level one when n q is large as the number of sparse 
grid points 7y(L, n q) is increasing exponentially with n q and L. In other words, L = 2 is preferred when 
SPDEs with many noises (large q) are considered. 

As discussed in the beginning of this section, the ratio of time cost for SCM and WCE is ry(L, ni/)/( N ^ n<? ). 
The cost of recursive multi-stage SCM with L = 2 is at most 1.8 times (1.6 predicted by the ratio above, 
q = 2 and n = 1) of that of recursive multi-stage WCE with N = 1. However, in this example, the 
accuracy of the recursive multi-stage SCM is one order of magnitude smaller than that of the recursive 
multi-stage WCE when N = 1 and L = 2. In Table 5, we present in bold the errors between 3.5 x 10 -5 
and 8.0 x 10 -5 . Among the four cases listed in the table, the most efficient, for the given accuracy above, 
is WCE with N = 2, which outperforms SCM with L = 3 and L = 2. WCE with N = 1 is less efficient 
than the other three cases. We also observed that when eri = (72 = 1, SCM with L = 2 is one order of 
magnitude smaller than WCE with N = 1 (results not presented here). 

For non-commutative noises in this example, we show that the error for WCE is A 2 + C A/n and the 
error for SCM is A. The numerical results suggest that SCM with L = 2 is competitive with WCE with 
N = 1 for both small and large magnitude of noises if n = 1. 

Table 6. Algorithm 2.2: recursive multi-stage Wiener chaos expansion for (4.7) at 
T = 1: ay = 0.5, (72 = 0.2, /3 = 0.1, e = 0.02. The parameters are M = 20, N = 2, and 
n = 10. The time step size 5t is A/10. 


A 

&\T) 

order 

-r.OO /rji\ 

<?2 ( T ) 

order 

CPU time (sec.) 

1.0e-l 

4.9310e-05 


2.6723e-05 


84.00 

5.0e-2 

1.4031e-05 

A 1 - 81 

7.3571e-06 

A 1.86 

160.50 

2.0e-2 

2.9085e-06 

A 1 - 71 

1.4171e-06 

A 1.80 

391.40 

1.0e-2 

9.8015e-07 

A 1 - 57 

4.4324e-07 

A 1.68 

749.40 

5.0e-3 

3.5978e-07 

A 1 - 45 

1.5082e-07 

A 1.56 

1557.60 

2.0e-3 

9.8910e-08 

A 1 - 41 

3.8369e-08 

A 1 ' 49 

3887.50 


With these three examples, we observe that the convergence order of the recursive multi-stage SCM in 
the second-order moments is one for commutative and non-commutative noises. We verified that our error 
estimate for WCE, A N + A 2 , is valid for commutative noises, see Examples 4.1 and 4.2; the numerical 
results for non-commutative noises, see Example 4.3, suggest the errors are of order A N + CA/n where 
C is a constant depending on the coefficients of the noises. 

For stochastic advection-diffusion-reaction equations, different formulations of stochastic products 
(Ito-Wick product for WCE, Stratonovich product for SCM) lead to different numerical performance. 
When the white noise is in the velocity, the Ito formulation will have stronger diffusion than that in 
the Stratonovich formulation in the resulting PDE. As stronger diffusion requires more resolution, the 
recursive multi-stage WCE with N = 1 may produce less accurate results than those by the recursive 
multi-stage SCM with L = 2 with the same PDE solver under the same resolution, as shown in the first 
and the third examples. 

To achieve convergence of approximations of second moments with first-order in time step A, we can 
use the recursive multi-stage SCM Algorithm 2.5 with L = 2, n = 1 and also the recursive multi-stage 
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WCE Algorithm 2.2 with N = 1, n = 1, as both can outperform each other in certain cases. For 
commutative noises, Algorithm 2.2 with N = 2 is preferable when the number of noises, q , is small and 
hence the number of WCE modes is small so that the computational cost would grow slowly. 

We also note that the errors of Algorithms 2.2 and 2.5 depend on the SPDE coefficients and integration 
time (cf. theoretical results of Section 3). For some SPDEs, the constants at powers of A in the errors 
can be very large and, to reach desired levels of accuracy, we need to use very small step size A or develop 
numerical algorithms further (e.g., higher-order or structure-preserving approximations, see such ideas for 
stochastic ordinary differential equations e.g. in [26]). Further, in practice, we need to aim at balancing 
the three parts (truncation of Wiener processes, functional truncation of WCE/SCM, and space-time 
discretizations of the deterministic PDEs appearing in the algorithms) of the errors of Algorithms 2.2 
and 2.5 for higher computational efficiency. 


5. Proofs 

5.1. Proof of Theorem 3.1. The idea of the proof is to first establish an estimate for the one-step 
(A = T ) error where the global error can readily derived from. We need the following two lemmas for 
the one-step errors. Introduce (cf. (2.6)) 

U N (t,x)= ^2 -/=< Pcc ( t , x )€ a - (5.1) 

|a|<N,aGj, y '°" 


Lemma 5.1. Let q = 1 in (2.1). Assume that 0 *,i, &*, c, i/i belong to Cf +1 (T>) and uo € H r (T>), 

where r > N + 1. Let u in (2.6) be the solution to (2.1) and un is from (5.1). For C\ < 8c, the following 
estimate holds 


E[|MA, •) - un(A, -)ll 2 ] < (C LrJ A) n +V^ a [ 


|A 


(C W A)W- 


N —1 


Sc 


(N + l)! 


L r J • 


8c-C! 


IKII 


flW 5 


where the constants Sc and Cc are from (2.3) and C i r j is from (3.5). 


Lemma 5.2. Under the assumptions of Lemma 5.1 and r > N + 2, we have 

PA 3 

E[||u N .n(A, •) - TtN(A, -)|| 2 ] < — 2 -C(N + 2,£)C(N, C)C N+2 e 2C ^ A+2C ^ A ||u 0 ||^, +a , 

n7r 

where Cc is from (2.3), C(N +2,£) is from (3.1), (7(1X1, £) is from (3.4) and CVj + 2 is from (3.5). 


Using Lemmas 5.1 and 5.2, we can establish the estimate of the global error stated in Theorem 3.1. 
Specifically, the one-step error is bounded by the sum of E[(u(A) — wn(A)) 2 ] and E[(un(A) — Ui\i,n(A)) 2 ], 
which are estimated in Lemmas 5.1 and 5.2. Then, the global error is estimated based on the recursion 
nature of Algorithm 2.2 as in the proof in [21, Theorem 2.4], which completes the proof of Theorem 3.1. 

Now we proceed to proving Lemmas 5.1 and 5.2. Let us denote by s k the ordered set (si, • • • , Sk) and 
for k > 1, denote ds k := ds\ ... dsk , and 


r(*0 


I = 

[ (■■■)ds k = 
J(k) 


ff 

ff 

Jo J si 


(• • •) ds i... dsk 




(■ ■ ■)dsk ■■■ ds 2 ds i, 


and F(A;s k ;x) = Ta- SI c M ■ ■ ■T S2 - Sl MT Sl u 0 (x), where M := M\. 

Proof of Lemma 5.1. It follows from (2.3) and the assumptions on the coefficients that (3.1) and (3.2) 
hold, cf. [11, Section 7.1.3]. Also, by the assumption that 0 ^ 1 ,iq belong to Cf +1 (T>), it can be readily 
checked that (3.3) holds. 
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By (2.6), (5.1) and orthogonality of (see (2.7)), we have 

E[|KA,.)-„ N (A,.)f]= ■£ E lly ° ( ^' )l12 - 

k>N \a\=k 

Similar to the proof of Proposition A.l in [21], we have 


E 

\at\=k 




p(k) 

/ \F(A;s k -,x)\ 2 ds k . 


Then by the Fubini theorem, 




|a|=/c 


(5.2) 


Assume that r > 0 is a integer. When r > 0 is not an integer, we use |_tJ instead. 

Denote X k = T Sk -s k -iM ■ ■ ■ T S2 - Sl MT Sl u 0 , Y k = MX k , k > 1 and also X = TA- Sk Y k . Then 
X k = 'T' Sk — Sk _ 1 \ k —\ and 1 k ~i = XlX k —i. 

By the dehnition of F, (3.1), (3.3) and (3.5), we have for r > k: 


||F(A;s fe ;-)|f < e 2CAA ~ sA \\Y k f H0 = e 2CAA ~ s ^ \\M lk X k \\ 


2 

H° 


H 1 


< C{Q,M)e 2Gc ^- Sk) \\X k \\ 

< Cl e 2CAA ~ s Hn-rll^ < • • •< C k e 2C ^ A \\u 0 \\ 2 Hk , 

where C{r — 1 ,M) is from (3.3) and C k is defined in (3.5). We then have 

r(*0 „ „ r( k ) 


> 9 n 1 

I \\F(A;s k ;-)l ds k <C k e 2C £ A \\u 0 \\ 2 Hk J ds k . 


(5.3) 


If r < k, by changing the integration order and applying (3.1), (3.3) and (3.2), we get 

dfc) . „ . /■(*) 9 r 

||W|| 2 ds k = / 

J(k) 


J y ||-f(A;s fe ; ■)|| 2 ds k = / 

< [ e 2C c( A ~s k ) ||y fc ||2 ds k = [ e 2Cc(A-s k ) \\ Mik X k \\ 

J(k) J{k) 

< C(0,M) [ e 2CAA ~ s ^ \\X k \\ 

J(k) 


All 2 ds k 


2 ds k 


< Sl'Ci 


\ H i ds 

0 2Cc(A — s k ) 

s 2 C £ (A-Sfc-i) ||y fc _ 1 || 2 ds k ~\ 


(5(0,M) f [ e 2UAA ~ s ^ \\X k \\ z Hl ds k ds 

J (k-l) J s k _i 


k-1 


l(k-l) 

where C\ is from (3.5). Repeating this procedure and using (5.3), we obtain 

r(fc) 


r(k) r />(r) 

j \\F(t;s k ;-)\\ ds k <6 r c ~ k d[- r J e 2CAA ~ s ^ ||y r || 2 ds r < 5 r c ~ k C k - r C r r e 2Cct ||u 0 ||^ r j ds r . (5.4) 





COMPARISON OF WIENER CHAOS AND STOCHASTIC COLLOCATION 


19 


By (5.1), (5.2), (5.3) and (5.4), and / ds k = —— , we conclude that, for r > N + 1 and C\ < Sc, 
E[|| U (A,.)-u n (A,.)|| 2 ] = E / ||F(A;s fe ;.)|r^ fc +^ / \\F (A; s fc ; •) || 2 ds k 


N<fc<r 1 


< E ^C k e 2C - A \\u 0 \\ 2 Hk + ^-Cfe 2C - A |M| 2 ^ E 


N<fc<r 


°Cr A (C' r A ) r " N ” 1 
L (N + 1)! J_ 


< (C r A) N + 1 e 2C ^[ 7 L_+ l ° rZ1J ' &C 1 "-" 2 


fc>r 

IM '% r n 


r! S c - G x 

Remark 5.3. Lemma 5.1 holds for r = oo if Coo < oo. Based on (5.3), we can prove that 
E[|| M (A, •) - m(A, -)H 2 ] < E T[ : - ck oo e2CcA IMI^ < (O.Af+V^E-E- IM*~ • 

fc>N ' ^ ’’ 

If r < oo, we need to require that C\ < Sc, i.e., (7(0, A4)C(1, £) < Sc- For example, C = A, 
A4i = \Dx, for which (7(0, Al)(7(l, £) = \ < Sc = 1. 


Proof of Lemma 5.2. It can be proved as in [21, p.447] that 


E[K(A,.)-UN,n(A,.)| 2 ] = E E E 

Z>n+1 k— 1 \ot\—k,i^—l 


¥&(A,-) 


(5.5) 


where i ^ a | is the index of last non-zero element of a and the last summation in the right hand side can 
be bounded by, see [21, (3.7)], 


E 


(A, x) ^ [V'-V * f s ^ 


< 


\cx\=k,i k =l 


/ (K-L) rs j+1 

E / F j (A;s k ;x)M l (s j )ds J 

A — 1 J S A — 1 


1=1 ^-i 


ds k , 


where ds k = dsi ■ ■ ■ dsj- 1 dsj+ 1 • • • dsk, sq := 0, s^+i := A, Mi(t) = / 0 * mi(s) ds and 

Fj{ A;s k -,x) = dF ^ s Sk ’ x) = TAs k M---Ts j+ 1 - Sj MCT Sj s j _ 1 ---T Sl uo(x) 


-TA-s k M ■ ■ ■ MCT Sj+1 - Sj ■ ■■T Sl u 0 (x) -. Ff +Ff. 


Then by the Fubini theorem and the Cauchy-Schwarz inequality, we have 


E 


iMAoir 


/ (fc-l) fc fSj+i 2 r s j+i 

E/ ||L(A; s fc ;-)|| dsj / M?( Sj ) d Sj ds k . 

■ i J Sa — 1 J Sa i 


|a|=fe,ig=J 

We claim (see its proof below) that 


7=1 Js o -1 


|F.,'(A;s fc ;-)ll < 2 max ||f/ II< 2CL 2 C(k, C)C{k + 2, C)e 2CcA ||u 0 ||^+2 . (5.6) 

l<j<k J 


Thus, by (5.6) we have 


E IIV?a( ^’ ,)l1 < 2kACk +2 C(k + 2, C)C(k, C)e 2CcA |K||^ +2 E M?(s)ds E ^ ds k . (5.7) 

\a\=k,i%=l 01 ' J 


r(k- 1) 
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Then by (5.5), (5.7) and Mi(t) = sin( 11 £ >7T t ) (by (2.12)), we obtain that 


d-i)” 


E[||u n (A, •) - M|\|,n(A, ■ 


* E 

Z>n+1 

2A 3 


A 2 N 0h/\ k 

e 2 C £ A ^ c k k+2 C(k + 2, C)C{k,C) \\u 0 \\ 2 Hk+2 


(l- 1) 2 7T 2 

N 


k =1 


(fc-l)! 


< 


< 


n7r z 

2A 3 

tm 2 


e 2C,A J- C k +2 C(k , £)C(fc + 2, £) ||uo||^ fc +2 kA 


k -1 


fc=l 


C n+ 2 C(N + 2,£)C(N,£)e 2CN + 2A+2CcA ||uo|| 


(*- 1 )! 

2 

H N + 2 * 


It remains to prove (5.6). Note that it is sufficient to estimate ||-F/|| due to the same structure of the 
two terms in Fj ( A ; s k ; x ). By the assumption that a % j and c belongs to C k + 3 ( V ), it can be readily 
checked that (3.4) holds with l < N + 1. Repeatedly using (3.1) and (3.3) gives 


\ F j\ 


MO 


= j|7A_ Sfc A4 • • • Ts j+1 - Sj M£T Sj -s j _ 1 • • • xto|| 

< e 2C ^ A ~ s ^ \\M---T Sj+ 1 - SJ M£T Sj - Sj _ 1 ---T S] 

< C(0,M)e 2Cc( ' A ~ Sk ' 1 || T Sk s k _ 1 ---T Sj+1 - Sj MCTs j - Sj _ 1 ---Ts 1 

< C ie 2C ^ A - Sk -^ \\M---T Sj+ 1 - Sj MCT Sj - Sj _ 1 ---T Sl u 0 \ 

< ■■■< c k k Z]e 2C ^ A ~ s ^ || MCTs.-s,! 


MO 


m 


h 1 


MO 


H k ~i 


\2 

I H k -i+ k 


< C k ~_]C(k - j,M)e 2C ^ A ~ s ^ || CTsj-s,-! ■ • -r si M 0 | 

< C k k :]C{k - j,M)C(k - J + l,£)e 2C ^ A -^ • • ■ T sl u 0 \\ 2 Hk _ j+3 

< C k k l£lC(k - j + l,£)e 2C ^ A -^ WTs.s^M • • -T S1 mo||^_ j+3 ■ 

where we have used (3.4) in the last but one line and the fact that C(k — j + 1, C) > 1. Similarly, we have 
||r Sj - 5j _ 1 M---r si Mo||^_, +3 < C{k-j + Z,£)Ci-\e 2C ^ ||mo||^ +2 . 

Thus, we arrive at (5.6). This end the proof of Lemma 5.2. □ 


5.2. Proof of Theorem 3.3. To prove Theorem 3.3, we need a probabilistic representation of the 
solution to (2.1). Let ({Rfe(s)} , 1 < k < d, Ff) be a system of one-dimensional standard Wiener processes 
on a complete probability space (ST 1 , J" 1 , Q) and independent of w(s) on the space (SI® SI 1 , F® J 71 , P®Q). 
Consider the following backward stochastic differential equation on (SI 1 , F 1 , Q), for 0 < s < t, 

d 

dX t ,x{ s ) = b(X t ^ x (s))ds + '^2a r (X t ,x{s))dB r (s), X t ,x(t) = x, (5.8) 

r— 1 


The symbol “d” means backward integral, see e.g. [19, 30] for treatment of backward stochastic integrals. 
The dx d matrix a(x) is defined by a(a;)a T (x) = 2a(x). Here a(x ) and b{x) are from (2.2). Consider the 
following backward stochastic differential equation on (SI ® SI 1 , F ® F 1 ,P <E> Q) for 0 < s < t, 


dYt,i, x (s) = c(X tjX (s))Y t ,i,x(s) ds + YjI=i Vr(X t:X (s))Y tAtX (s)\ dw r , Y t ,i, x (t) = 1. (5.9) 


Here c(x) and v r (x) are from (2.2). When it 0 (;r) €= C 2 (V) and a(x ), b(x), c(x ), v r (x ) S C°(2?) and 
oy r = 0, the solution to (2.1)-(2.2) can be represented by, see e.g. [19], 


u(t 7 x) = EQ[w o (A!t,x(0))exp( 


E 

r=l ' 


v r (X t , x (s)) dw r (s) 


c(X M (s))ds)], 


(5.10) 


where c(x) = c(x) — \ Er=i v r( x )- 
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Here we first establish the one-step error (3.9) and then the global error (3.10) We follow the recipe 
of the proofs in [16, Theorem 3.1] and [5, Theorem 4.4] where n = 1 and K > 1. 

We need the following mean-square convergence rate for the one-step error. 

Proposition 5.4 (Mean-square convergence). Assume that ai^ r = 0 and that the initial condition u 0 and 
the coefficients in (2.2) are in Let u(t,x) be the solution to (2.1) and u A , n {t,x) the solution to 

(2.16). Then for any e > 0, 

E[|u(A, x) — u A >n (A,a:)| 2 ] < C exp((7A)(A 3 + A 2 )n _1+e , (5.11) 


where the constant C > 0 is independent of n. 

Proof. The solution to (2.16) using the spectral truncation of Brownian motion 4 A,n * from (2.13) can 
be represented by, see e.g. [5, 16], 

u A ,n(A,a) = E Q [Mo(X Aja; (0))exp(X;^ = i/ 0 A ^r(^A, :E (5))rf'«;.r A ’ n) (5) + / 0 A c(X Ai3 .(s))(is)]. 

(5.12) 

Using e x — e v = e ex+< J~ e ' )v (x — y), 0 < 9 < 1, boundedness of c{x) and uq{x) and the Cauchy-Schwarz 
inequality (twice), we have for some C > 0: 


JwA,n(A, x) - u(A,x)| ] 


(5.13) 


/ pA 9 pA 

E[( E Q [uo{X AtX (0))exp(J c(Xa,x{s)) ds )) exp(y^ J iy r (X A>x (s))[9 dw^ A ’ n \s) + (1 - 9) dw r (s)]) 


< 


9 p A 

(V / t / r(^A,a(s))[(iui^ A - r ' ) (s) - dw r (s)])] 

r=l Jo 

9 r A 

£rdexp( 


< 


C exp(CA)E[ |e q [exp(y^ J v r (X AtX (s))[0 dwl A ' n \s) + (1 - 0) dw r {s)]) 

9 pA 

V / v r (X A , x (s))[dwl A ' n) (s) - dw r (s)] 

r=l J ° 

9 pA 

C exp(CA) (E[Eq[ exp(E' / 4:V r (X AtX (s))[9 dw( A,n \s) + (1 - 0) dwRs)])]]) 1 

r=l J ° 

x(E[E QlC^ j A v r {X A ^{s))[dw^ n \s) - d Wr ( S )]) 4 ]]) 1/2 . 


Recall that E[-] = Ep[-] is the expectation with respect to P only. Hence, we need to estimate Ij = 
(E[E Q [(E?=i/o A ^(^A,,(s))[d^ A ' n) ( S )- dw r (s)}) 4 }]) 1/2 and 


I 2 = 


(E[(Eq [exp(^ J 4v r (X A , x (s))[6 dw( A ’ n \s) + (1 - (9) dw r (s)])]]) 1/2 . 


We first estimate R. Due to the independence of Bk and w r , and according to [28] and (2.5), we have 
pA pA pA 

/ v r {X A ^(s))dw r {s) = / Vr{X A , x {s))odw r {s) = Y j ti r ,i / v r {XA, x {s))m r ^{s)ds. 

Jo Jo i=0 Jo 
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Thus by the Fubini theorem, (2.5) and (2.13), we can represent I| as 


U — (EglM 

= (Eq[E[ 


9 /* Ai /* Ai 

v r {X A iX (s))du;(. A ’ n) (s) - / i/ r (i Ail (s))orfw r (s)] 
_i Jo Jo 


r— 1 1/0 
q 00 


]]) 


1/2 


y ^ p LA 

^ ^ ^ £,r,i I ^ 

;_, 1 ^0 


r=l i—n+1 
<? 00 r A 


]]) 1/2 


< 


H _ _ /* A\ 

( 3E q[(X^ W ^(^A^(5))m r ,i(5)ds) 2 ) 2 ]) 1/2 , 

^-1 -n I 1 0 


r—1 2—n+1 


where we have used twice the fact that X A ,i are independent of w r and 1 • Then by standard 
estimates of L 2 -projection error (cf. [7, (5.1.10)]), we have for 0 < £ < 1, 

oo r A 2 


C-VJ p / \ 

y" (/ ^ r (i'A,x(s))m T . ii (s) ds) 2 < C'A 1_e n _1+E is r (X AtX (-)) 
, do 


2 =n+l 


-,2,[0,A] 


(5.14) 


where the Slobodeckij semi-norm \f\ 6p j 0 A j is defined by (/ Q A / A dxdy ) 1 / p and the constant 


A 1 £ appears due to the length of domain, see e.g. [7, Chapter 5.4]. Thus, we obtain 


Ii < CA 1 - £ n- 1 + E (y]EQ[ MX A A')) 


, ]) 1/2 , 0<e<l. 

A=i,2,[0,A] Jy 


(5.15) 


By (5.8) and the Ito formula, we have 

ps P 

Xa,x{ S ) - XA,x(ai) = / b(X A ,x( S 2)) dS2 + (-^A,a:(si))[JBfe(s) - -Bfc(si)] + R(si,s), 

d«i fc—i 


where Eq[|.R(si, s) | 2Z ] < C|si — s| 21 (l > 1) when b{x) and ak(x) belong to C 2 (X>). By the Lipschitz 
continuity of ly, the definition of the Slobodeckij semi-norm, it is not difficult to show that 

4 


1 21 


E 


■Q\ 


Xa,x(')) 


i=i,2,[0,A] 


] < C(A 4+2e + A 2+2e ). 


Thus, by (5.15) and (5.16), we have 


Ii < C(A d + A 2 )n 


2\„-l+e 


(5.16) 


(5.17) 


Now we estimate I 2 . Using the following facts (see e.g. [16, Lemma 2.5]) 
<1 r A 1 ,A 


H f‘ LA Sf PlA 

E[exp(V' / 4:is r (X A .x(s))dw r )\ = exp(Y] 8 / A(^a,x(s)) ds), 

r=l J 0 r—1 do 

9 „a q r a 

E[exp(^ J 4 v r {XA,x(s))dwl A ’ n) {s)} < 4exp(^8 J is?(X A , x (s)) ds), 


we have I 2 < 4exp(CA). From here, (5.17) and (5.13), we reach (5.11). □ 

Now we are ready to prove Theorem 3.3, i.e., the convergence in the second moments. For simplicity 
of notation, we consider q = 1 while the case q > 1 can be proved similarly. Denote 

rt 

m > n. 

2=1 j=n+l 


n m „t 

UA, n , m ,s{t,x,y) =: uo(W,:z(0))exp(V'v M y i + 6> V] v x ,jVj + / c(X M (s)) ds), 

_, , do 


where = / 0 * (X tiX (s))mi(s) ds for * < m (X t ,x(s) is the solution to (5.8)) andy = (yi ,..., y n ,y n+1 ,. 

Let us write UA, n . m ,e(t,x,Z) = ^Q[U A ,n,m,e{t,x, S)], where S = (£ 1 ,... , £ n >£n+i> • • • ,£m)- With this nota¬ 
tion, we have 

RA,m(t,^) — fiA,n,m,l (t, X, Zl) , WA,n(t,Ic) = UA,n,m,o(^; ^)- 


■ • i 2/m)- 
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For m > n, by the first-order Taylor expansion, we have 
| E [wA,m( A ^) - «A,n( A 5 a ; )]| 


= |2 E 

i,j— n+1 

m 

+2 E 

i,j=n+l 


1 


($i,j + 1) JO 

1 t 


9(1 - 6»)E[uA,n,m,fl(A, x , S)Eq[C/A, n ,m,e(A, x, E)u lti (t, x)vij(t, x)]&£ 3 -] d-0 


(8iJ + 1) J Q 


9(1 - 0)E[Eq [f7 A ,n,m,e(A, x, E)i x)]Eq [U A , n , m ,e(t, x, S)u ltj (t, x)]&£j] d9\ 


< 2 


[ (1 - 0)0E[uA, n .m,e(A,x,E;)EQ[17A,n, m ,e(A,x,E:) [ Y^ ^i,i(A,x)^ ] }}dd 

J ° \i=n+l ' 

J o (1-0)0E[^£] EQ[£/A,n 1 m,fl(A,®,HK i (A,s)]^ ]d0 


(5.18) 


where Si tJ = 1 if i = j and 0 otherwise and we have used the facts that * > n, are independent of 
u n (t,x) and E[£j] = 0. 

By the Cauchy-Schwarz inequality (twice), we have for the first term in (5.18): 

J (1 - d)6»E[M A , n , m , e (A,x,S)EQ[17 A:n , m: e(A,x,S) ( ^ zq^A^&J ]\dd 


i=n+l 


< C(E[E Q [ ]T ^1,*(A, x)£i j ]]) 1/4 . 


(5.19) 


yi=n+l 
\~2 


Here we also used that E[u A n , m g( A, x, S)],E[(EQ[{7 A . n m,e(A, x, S)]) 2 ] < C, which can be readily checked 
in the same way as in the proof of Proposition 5.4. 

By the Taylor expansion for b r A.n.m, 6 »(A, x, y), we have 

m „i 

£/A,n,m. 0 (A,x,y) = t/ A ,n,m.o(A,x,y) + Y] u lti (A,x) / (1 - 0i)0i0£/A,n,m,e0 1 (A,X, y) dOiyt- 

i—n+l 

Then by the Cauchy-Schwarz inequality (several times) and the fact that £j, i > n are independent of 
f^A,n.m,o(tj x, 5), we have for the second term in (5.18): 

2 

^ .. .. dO 

\i= n+l 


< 4 


J Q (1-0)0E[^ EQ[C/A,n,m, 0 (A,X,S)^ >i (A,x)]^j ], 

pi m 

/ (1-0)0 ]T E[(E Q [C/ A ,n,m,o((A,x,S)^ >i (A,x)]) 2 ]d0 

i= n+l 

f (1-9)9 3 E[(e q [[ (l-0 1 )0 1 J7 A .n.m.00 1 (A,x,H)d0 1 ( V ^(A,^)) 2 ^ ]d0 

• / 0 V 170 i=n+l / 

m 

< E[E Q [C/ 2 , n , mi0 ((A,x,5)]] £ Eq[i/^(A,x)] 

z=n+l 

<■ 1 /*1 m 

(1 — 0)0 3 (E[Eq[ / (l-0 1 )0 1 f/l. n!m ^ 1 (A,x,5)d0 1 ]) 1 / 2 d0] (E[E Q [( ^ p m ((A,x )&) 8 ]]) 1/2 
) Jo - „ , 1 


i=n+l 


< E Q [^(A,x)] + C(E[E Q [( ]T MA,x)&) 8 ]]) 1/2 . 


(5.20) 


z=n+l 


i=n+l 
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Here we used that E[Eq[C/^ n m 0 ((A, x, £)]], E[Eq[{ 7^ n m (A, x, S)]] < C, which can be readily checked 
in the same way as in the proof of Proposition 5.4. 

By (5.18), (5.19) and (5.20), we have 

| E [“A,m( A ’ X ) -«A,n( A ^)]| ( 5 - 21 ) 

mm m 

< c ^QWli(A,x))+CE[E Q [( E ^ 1 , i (A,x)^) 8 ]]) 1/4 + C(E[E Q [( E id,i(A,a:)£») 8 ]]) 1/2 . 

i=n+l i=n+l 2 =n+l 

Similar to the proof of (5.15), we have 


m / m \ 4 

E[Eq[(£ ^.((A,®)^) 8 ]] < CEq[( J2 ^m( A ^) }<CA 4 ^^ 1 -^E Q {\v 1 (X AtX (-)) 

\i= n+1 / 


i=n+l 


^. 2 . 10 , A] 


Similar to the proof of (5.16), we can estimate Eq[ ^(Xa^Q) 

Eq 


i^,2,[0,A] 


1^2,2,[0,A] 

] < C(A 8+4e + CA 4+4e ), 


as follows: 


and thus 


Similarly, we have 


^i( a a,,(-)) 

m 

E[Eq[( E (A,x)&) 8 ]] < C(A 12 + CA 8 ). 


i— n+1 


E[Eq[( E M a .*)&) 2 ]] = E EQKi(A,x)]<C7(A 3 + CA 2 ). 


(5.22) 


(5.23) 


2 =n+l 


i=n+l 


By (5.21), (5.22) and (5.23), we have 

| E [“A,m(A,a:) - UA, n (A,a;)]| < C'exp(CA)(A 6 + A 2 )rT 1+E . (5.24) 

By the triangle inequality and the Cauchy Schwarz inequality, we obtain 

|E[u 2 (A,x) - UA, n (A,a:)]| < | e [m 2 (A,x) - «A, m (A,®)]| + |E[hi m (A,a;) - ui n (A,x)]| , 

< C(E[|u(A, x) - UA,m(A, x)| 2 ]) 1/2 + |E[ui m (A, x) - ui >n (A, x)] | . 

The one-step error (3.9) then follows from (5.24), Proposition 5.4 and taking m to +oo. The global error 
(3.10) is estimated from the recursion nature of Algorithm 2.5 as in the proof in [21, Theorem 2.4]. □ 


5.3. Proof of Theorem 3.4. For any n-dimensional function <p(yi, ... ,y n ), we denote 

L ,p{yi ■ ■ ■ • 1 »"> ■ exp (" \ g ) dy - 

Introduce the integrals 

EV = / ^C^ 1 ’' ‘' ,Vk '' • ■ ’^ n ) ex P dykl fc==1 >"-’ n ’ ( 5 - 25 ) 

(k) 

and their approximations Qn ; by the corresponding one-dimensional Gauss-Hermite quadratures with n 
nodes. Also, let — QiP-v By the definition of Smolyak sparse grid and using the recipe from 

the proof of Lemma 3.4 in [27], we obtain 

n 

In<P - A(L, n)ip = E S(L, n) ®£ = m l[ k \ + (J< 1} - Q< 1} ) ®U 
1=2 


(5.26) 
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where 


S(M)= E -fc=i -4 

iiH-Mz—i"Mz 


{-JM? ® ( 4 ° - q 1 °). 


(5.27) 


Denote by D a the multivariate derivatives with respect to y. According to the proof of Proposition 
3.1 in [38], we have 


< 


5(L,o®^ +1 /rv 

(3c/2)# Gi - 1+1 


(5.28) 


E 


uH-Mz—L+Z—1 




® mG F i _ 1 gE^ 2 “v( y) 


^GGz_i 


exp[- y f-f-EII n 


»2 


/c—Z —|— 1 1 n£.Gi ^ 


where the multi-index cxj = (ii — 1,..., i/_i — 1, ip 0,..., 0) with the m-th element a™, the sets i = 
Ji_i(oi) = {m : a™ = 0, m = 1,..., Z — 1} and G;_i = Gi-i(ai) = {m : a™ >0, m = 1,..., l — 1}, 
the symbols #A;-i and #Gj_i stand for the number of elements in the corresponding sets. Here c > 

0, 0 < /3 < 1 are only related to the Gauss-Hermite quadrature Q , and are independent of the number 
of nodes in the Gauss-Hermite quadrature, see e.g. [24, Theorem 2]. 

Proof of of Theorem 3.f. Setting tp(y±, • • • , y n ) = u\ n (t, x,y i, • • • , y„), we then have that A(L, n )ip is 
an approximation of the second moment of the solution obtained by the sparse grid collocation methods. 

Recall from (5.12) that u A , n (t,x,y) = EQ[U A , n (t,x,y)] where U A ^ n {t,x, y) = U A . n , mfi (t, x, y). 

Now we estimate D 2ai [u\ n (A, x, yi, ■ ■ ■ ,y„)]. To this end, we need to first estimate D /3i [uA, n (A, x, y\, ■ ■ ■ , y n )\, 
where fa < 2cq. By (5.14), we have for 0 < e < 1, 


x ) — G(Amax(fc — 1, l)) e 1 vi{X A , x (-)) 
we have, by the Cauchy-Schwarz inequality, 


2 

i=S,2,[0,A] 


|D /3! UA,n(A,a;,y)| = 


EQ[G A .n(A,a;,y) ]J(pi ifc (A, x))^ 


k=l 


(5.29) 


< ( E Q[C / A.n(A,x,y)]) 1/2 (E Q [J|(i/i )fc (A,x)) 2ftfc ]) 1/2 


k=1 


k =2 


< (GA 1 - £ )lftl/ 2 [](fc - l) (£ - 1)ftfc/2 (E 0 [Gi in (A,x,y)]) 1 / 2 (E Q [ ^(X A ,,(-)) rn J) 1/2 - 


-,2,[0,A] 


By the chain rule for multivariate functions, we have 


7? 2 “‘ [«i,n(A, x, y)] = E ( 2a ') ! 

A+7l=2ai 


D^ l u A . n (A,x,y) D 71 UA, n (A,x,y) 


A! 7i ! 

and thus by (5.29) and the fact that ^ft+ 7i = 2 aj = 2 2 l ai l _1 , we have 

i 

\D 2a ‘[ul^A,x,y)}\ < 2 2 ^-\CA 1 -^E Q [Ul n (A,x,y)} (k - l) (e - 1)a ? 


k—2 


x max ((Eq[ ui(X AtX (-)) 

A+ 7 l=2aj 


2|A| 


iA,2,[0,A] 


]) 1/2 (Eq[^(4x(-)) 


4|:7i I 
iyS,2,[0,A] 


]) 1/2 ). 
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Similar to (5.22), we have Eq[ vi(Xa,x(-)) 


2|ft| 

i=i,2,[0,A] 


] < C{ Alftl( 2+e ) + and 


| D 2a ‘ [{&_„(A, x, y)] | < C(A 3 '“'I + A 2 '“‘l) ]J (* - l) (e - 1)af E Q [t/i n (A, x, y)]. (5.30) 

k—2 

Then by (5.28) and (5.30), we obtain 

s(u)®(u+i4 n V 

l 

< C( A 3L + A 2L )(1 + (3c/2) iAi ) / 9- (LAi)/2 E[E Q [C7i, n (A,x,y)]] E []> “ l) (s ~ 1)a ‘ 

i,\ —L—|— Z—1 /c—2 

< C(A 3L + A 2L )(1 + (3c/2) iAi )^- (LAi)/2 e 1 - L (/ - l) u_1 (5.31) 

with the constant C > 0 which does not depend on n, e, L, c, /?, and l. In the last line we used the fact 
that E[Eq [t^A.n(^> m ? y)]] is bounded and that 


e n (fc- i ) (e - i) “?=a - i ) e - 1 


ii -(-• • —L-|-Z—1 k —2 


E IK*- 1 )' 

—L+Z—1 k —2 


e—l)(*fc — 1) 


< (/ - l) e " 1 (E(fc ^ If 1 ) 1 ” 1 <(l- l) e-1 (e _1 (i - 1) £ ) L_1 = e 1_L (^ - I) 16 " 1 . 

k—2 

Then by (5.26) and (5.31), we have 

n 

\I n <p-A(L,n)ip\ < C(A 3L + A 2L )(l + (3c/2) iAn )/3" (LAn)/2 e 1 - L E(^ 

1=2 

+ |(A (1) - Q^) 

< C(A 3L + A 2L )(1 + (3c/2) LAn )/3-( LAn )/ 2 £- L L- 1 n Le , 

where the term in the second line is estimated by the classical error estimate for the Gauss-Hermite 
quadrature Q , see e.g. [24], and the estimation of derivatives (5.30). 

The global error is estimated from the recursion nature of Algorithm 2.5 as in the proof in [21, Theorem 
2.4], □ 
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